################################################################################
#                                                                              #
#        Replication script for analysis contained in:                         #
#        BRIEF OF AMICUS CURIAE FAIR TRIAL ANALYSIS, LLC IN SUPPORT            #
#        OF WARREN KING’S PETITION FOR WRIT OF HABEAS CORPUS                   #
#        Filed Feb. 5, 2025 in Superior Court of Butts County, Georgia         #
#        Case No. 2024-SU-HC-0010                                              #
#        Judge: Hon. Robert L. Mack                                            #
#                                                                              #
################################################################################


################## Set-up machine and load datasets ############################

# install the following packages: sate, RCPA3
# download king_survey_data.csv and cloud_respondent_info.csv to working directory

rm(list=ls())  # starting with clean environment
library(sate)  # make sure you're using version 2.3.0
library(RCPA3)
# set working directory to location of downloaded datasets on your machine
setwd("C:/Users/Barry Edwards/Documents/warren king case") # edit for your machine

#### Load datasets
surveydata <- read.csv("king_survey_data.csv")
respondentdata <- read.csv("cloud_respondent_info.csv")


########################## Set-up for analysis #################################

#### Identify jury-qualified respondents

# qualifying the jurors
# Q10 asks about standard qualifications
# Q16 asks about excluded occupations
# Q18 asks about following death penalty instructions
isQualified <- (surveydata$Q10=="Yes" & surveydata$Q16=="No" 
                & (surveydata$Q18!="It should never be used. I would never support putting someone to death regardless of what they did." &
                   surveydata$Q18!="Murderers should always be put to death regardless of how sorry they are now or how rough their life has been."))

## calculate sampling weights for analysis
# encode Cloud Research respondent data in appropriate format
respondentdata <- encode.cloud.respondent.variables(dataset = respondentdata)

# define the demographics of the target population (i.e. nationally representative)
target.demog <- target.population.demographics("GA")
# not using GA level data, using Butts Co. data instead
target.demog$ba_or_more <- .247 
target.demog$woman <- .463 
target.demog$hispanic <- .036
target.demog$black <- .276 
target.demog$hhincome_over50k <- .607 
target.demog$age35plus <- 14477 / 20724 

respondentdata <- weights_for_population(respondentdata = respondentdata, targetdata = target.demog)
####  process the respondentdata for weighting before merging it with survey data
surveydata <- merge(surveydata, respondentdata, by = c("ParticipantId"), all.x=TRUE)
# there may be some survey responses that don't match to a respondent in merged dataset
# to be safe, code those survey responses with mean weight  
surveydata$weights[is.na(surveydata$weights)] <- mean(surveydata$weights, na.rm=T)
# adjusting for losing observations from jury qualifications, make sure mean = 1
surveydata$weights.qual <- surveydata$weights * (1 / mean(surveydata$weights[isQualified]))
# weights.qual is the variable used for weighting observations 


########################## Analyze data ####################################

## check representative of raw, weighted, and qualified & weighted samples
## results reported as Table 3

demog.variables <- c("ba_or_more", "woman", "hispanic", "black", "hhincome_over50k", "age35plus")
# Butts County, Georgia is reference (target.demog)
# weighted sample
fun <- function(x, w) wtd.mean(x, w=surveydata$weights)
apply(X = surveydata[, demog.variables], FUN = fun, MARGIN = 2)
# unweighted sample
fun <- function(x) wtd.mean(x)
apply(X = surveydata[, demog.variables], FUN = fun, MARGIN = 2)
# qualified & weighted sample
fun <- function(x, w) wtd.mean(x, w=surveydata$weights[isQualified])
apply(X = surveydata[isQualified, demog.variables], FUN = fun, MARGIN = 2)


##################### Analyze and Compare Sentencing Preferences ###############

# Code respondent votes in actual and hypothetical conditions
# Q26 is vote in initial condition
# Q31 is vote in cross-over condition
surveydata$vote_in_actual <- surveydata$vote_in_hypo <- NA
surveydata$vote_in_actual[surveydata$initial_assignment=="Actual"] <- surveydata$Q26[surveydata$initial_assignment=="Actual"]
surveydata$vote_in_actual[surveydata$initial_assignment=="Hypothetical"] <- surveydata$Q31[surveydata$initial_assignment=="Hypothetical"]
surveydata$vote_in_hypo[surveydata$initial_assignment=="Actual"] <- surveydata$Q31[surveydata$initial_assignment=="Actual"]
surveydata$vote_in_hypo[surveydata$initial_assignment=="Hypothetical"] <- surveydata$Q26[surveydata$initial_assignment=="Hypothetical"]

# find P(g|actual) and P(g|hypo)
freqC(surveydata$vote_in_actual[isQualified], w=surveydata$weights.qual[isQualified], plot=F)
freqC(surveydata$vote_in_hypo[isQualified], w=surveydata$weights.qual[isQualified], plot=F)

# Compare juror stats, reported in Table 4
compare.juror.stats(pg_actual = .2705, n_actual = 786, pg_hypo = .1981, n_hypo = 786)

# Table 5: Comparison of votes in actual and hypothetical trial conditions
# note that rows and column ordered differently in brief for clarity
crosstabC(dv=surveydata$vote_in_hypo[isQualified], iv=surveydata$vote_in_actual[isQualified], 
          w=surveydata$weights.qual[isQualified], ivlabs = c("Death", "Life"), dvlabs = c("Death", "Life"), 
          digits=1, plot = F)


##########################  Set-up for Automated Text Analysis  ################

## output a .csv file of qualified respondent comments for automated text analysis
write.csv(x=data.frame(weight=c(surveydata$weights.qual[isQualified], surveydata$weights.qual[isQualified]), 
                       comment=c(surveydata$Q28[isQualified], surveydata$Q34[isQualified])),
          file = "king_survey_respondent_comments.csv")
## weighted content analysis of comments performed with AI, not R. 


########################## Analyze and Compare Verdict Probabilities ###########

compare.jury.stats(pg_actual = .2705, n_actual = 786, 
                   pg_hypo = .1981, n_hypo = 786, jury_n = 12,
                   pstrikes=15, dstrikes=15)
# Per OCGA 15-12-165, both sides 15 strikes in death penalty case

## create a graphic illustrated change in verdict probabilities
## this is Figure 2 in brief
png(file="effect on verdict probabilities.png", width=6.5, height=3.5, 
    units="in", type="cairo", pointsize=11, res=300, antialias="default")
graphics::par(mar=c(2.5,2.5,1.5,.5), mfrow=c(1,2), omi=base::c(0.1,0.1,0.1,0.1), family="serif")
basic.plot.grid(xlab="",ylab="", main="")
mtext(side=1, text = "Jurors' Support for Death Sentence", cex = .8, line = 1.5)
mtext(side=2, text = "Probability of Death Verdict", cex = .8, line = 1.5)
mtext(side=3, text = "(a) Effect with Neutral Summaries", line = .5, font = 2)
lines(x = seq(0,1,.01), y = as.jury.point(sample_pg = seq(0,1,.01), jury_n = 12, pstrikes = 15, dstrikes = 15))
sate:::plot.ellipse(pg=.2704, n=786, point.col="gray25", pstrikes=15, dstrikes=15)
sate:::plot.ellipse(pg=.1981, n=786, point.col="white", pstrikes=15, dstrikes=15)
graphics::legend(x=.50, y=.2, legend=base::c("Actual Trial", "Hypothetical Trial", "Confidence Interval"),
                 col=base::c("gray25", "black", "gray25"),
                 pch=base::c(16, 1, -1), lty=c(-1, -1, 3), 
                 cex=.7, pt.cex=base::c(1.2, 1.2, -1),
                 box.col="#FFFFFF60", bg="#FFFFFF60")
constant_shift <- .35
basic.plot.grid(xlab="", ylab="", main="")
mtext(side=1, text = "Jurors' Support for Death Sentence", cex = .8, line = 1.5)
mtext(side=2, text = "Probability of Death Verdict", cex = .8, line = 1.5)
mtext(side=3, text = "(b) Potential Effect with Races Identified", line = .5, font = 2)
lines(x = seq(0,1,.01), y = as.jury.point(sample_pg = seq(0,1,.01), jury_n = 12, pstrikes = 15, dstrikes = 15))
sate:::plot.ellipse(pg=.2704+constant_shift, n=786, point.col="gray25", pstrikes=15, dstrikes=15)
sate:::plot.ellipse(pg=.1981+constant_shift, n=786, point.col="white", pstrikes=15, dstrikes=15)
graphics::legend(x=.50, y=.2, legend=base::c("Actual Trial", "Hypothetical Trial", "Confidence Interval"),
                 col=base::c("gray25", "black", "gray25"),
                 pch=base::c(16, 1, -1), lty=c(-1, -1, 3), 
                 cex=.7, pt.cex=base::c(1.2, 1.2, -1),
                 box.col="#FFFFFF60", bg="#FFFFFF60")
dev.off()

## analyze consequences of likely racial dynamic (as constant shift)
constant_shift <- .35
compare.jury.stats(pg_actual = .2704+constant_shift, n_actual = 786, 
                   pg_hypo = .1981+constant_shift, n_hypo = 786, jury_n = 12,
                   pstrikes=15, dstrikes=15)




